//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// dgrid_rt.cpp - Do GRID Roche Tomography: Programa para mapear a resolucao espectral em 
// arquivos de dados e preparar imagens das posicoes de cada pixel em cada espectro em cada
// fase. 
//
// Formato de utilizacao: dgrid_rt dados.dat spec.lis output.dgrt
// dados.dat - Arquivo com parametros orbitais padrao (veja roch_tom.cpp)
// spec.lis  - Lista com o nome dos espectros no formato nome_fase, onde fase é a fase orbital
//             do espectro.
// output.dgrt - Arquivo de saida com o formato que deve ser lido pelo programa de tomografia.
//
// Tiago Ribeiro - 2007 LinuxC++
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>
//#include <gsl/gsl_errno.h>
//#include <gsl/gsl_math.h>
//#include <gsl/gsl_deriv.h>
#include "light_curve.h"
#include <gsl/gsl_blas.h>
#include "convolve.hpp"
//#include <gsl/gsl_multimin.h>
//#include "light_curve.h"
//#include "light_curve.cpp"
//#include "roche.cpp"
//#include "lib_zabot/dataFileReader.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"

using namespace std;

// Minhas estruturas
template <> void Roche::NewI <vector<double> >(vector<double> &newi)
{
  if(newi.size() != star.theta.size()) aborta("Novo mapa de tamanho diferente do esperado!");

  for(int i = 0; i < newi.size(); i++){
    star.I[i]=newi[i];
  }

}

struct DataParKeeper {
  //vector<vector<double> > LineProfVel;
  //vector<vector<double> > LineProfFlux;
  //vector<vector<double> > LineProfFerr;
  vector<Histograma> dados;
  vector<double> Fase;
  vector<vector<double> > in_par;
  vector<vector<int> > pixmap;
  surface ImgDef;
  surface ImgSig;
  double PFac;
  double Qui2;
  double S;
  double dX2;
  double dS;
  double GX2;
  double GS;
  Roche *estrela;
};

// Minhas funcoes
double Radii2Ffac(const vector<vector<double> >&,int,int);
int definicoes_obscurecimento(const int law);
void show_par(ostream& out, vector<vector<double> > &par);
vector<vector<double> > org_par(vector<vector<double> > &in_par); //organiza os

int main(int argc, char *argv[])
{
  // Testar argumentos de entrada...

	if(argc < 3 ){
    cout << "dgrid_rt.cpp - Do GRID Roche Tomography: Programa para mapear a resolucao\n" 
	 << "espectral em arquivos de dados e preparar imagens das posicoes de cada pixel\n" 
	 << "em cada espectro em cada fase.\n\n" 
	 << "Formato de utilizacao: dgrid_rt spec.lis dados.dat output.dgrt              \n"
	 << "dados.dat - Arquivo com parametros orbitais (rode com '-h' para detalhes)	 \n"
	 << "spec.lis  - Lista com o nome dos espectros no formato nome_fase, onde fase  \n"
	 << "é a fase orbital do espectro.						 \n"
	 << "output.dgrt - Arquivo de saida com o formato que deve ser lido pelo programa\n" 
	 << "de tomografia.							       \n\n"
	 << "Tiago Ribeiro - 2007 LinuxC++                                               \n";

		if(argc < 2) return 0;
		if(!strcasecmp(argv[1],"-h")){
		   
		   cout << "\nFormato do arquivo de parametros:\n"
		   << "\tSIZE Xaim FWHM S/R ABS_FLAG texp/Porb\n"
		   << "\ti M1 M2 P vmax lv\n"
			<< "\tf2 grav_coef limb_law nspot Tpole R2[Rsun] fd\n"
			<< "\t(coeficientes de obscurecimento de borda)\n\n"
			<< "\n\t -> FWHM = largura utilizada para gerar mapa default (~0.1).\n" 
			<< "\n\t -> ABS_FLAG = -1 p/ emissao e 1 para absorcao.\n" 
			<< "\n\t -> texp/Porb = considera blurring instrumental. 0 para nao fazer nada e [0:1) para considerar.\n"
			<< "\n\t Note que 1 significa que texp = Porb!"
			<< "\n\t -> Esta versao le o Raio da secundaria em Rsol e converte para ffac.\n";
		}		   

    exit(1);
  }

  // Le os perfis de linha e armazena em estrutura do tipo DataParKeeper.

  DataParKeeper InFiles;

  DataFileReader lis_in(argv[1]); //lis_in eh lista de arquivos!

  cout << "#\n#Lendo arquivo de dados de " << argv[1] << "...\n#" << endl;

  int LineFlag = 0;
  if(lis_in.lineNumber() > 0)
    LineFlag = lis_in.lineNumber();

  while(!++lis_in){ // le lista ate o fim
    string filename = lis_in.readColumn<string>(1);
    DataFileReader lpf(filename.c_str()); // abre arquivo da lista!

    //dar um geito de pegar fase do nome do arquivo e...

    string::size_type las = filename.find_last_of("_");
    string Sfase (filename, las+1);
    double fase = atof(Sfase.c_str());

    cout << "# " << las << " - " << Sfase << " = " << fase << endl;

    InFiles.Fase.push_back(fase);
    //InFiles.LineProfVel.push_back(vector<double>(0));
    //InFiles.LineProfFlux.push_back(vector<double> (0));
    //InFiles.LineProfFerr.push_back(vector<double> (0));

    vector<vector<double> > rdata(2);

    while(!++lpf){ // le arquivo aberto ate o fim
      if(lpf.numberOfColumns()!= 3){ 
	ostringstream erro; 
	erro << "Numero de parametros errado\n\tArquivo: " 
	     << filename << "\n\tLinha  : " << lpf.lineNumber();
	aborta(erro.str());
      }
      int VL = lis_in.lineNumber()-LineFlag;
      double datain[2];
      //lpf >> datain[0] >> datain[1] >> datain[2];
	  //cout << "# Making apropriate conversions...\n"; 
      datain[0]=lpf.readColumn<double>(1) * 1e5;
      datain[1]=lpf.readColumn<double>(2);
		//cout << datain[0] << ' ' << datain[1] << endl;
      //lpf.seteol();
      //guarda na estrutura...
      rdata[0].push_back(datain[0]);
      rdata[1].push_back(datain[1]);
      //.push_back(datain[2]);
      //++lpf;
    }
    Histograma rhist(rdata);
    InFiles.dados.push_back(rhist);
    cout << "#Arquivo " << filename << " pronto..." << endl;
    lpf.close();
    //++lis_in;
  }
  lis_in.close();
  int Nsize = 0;
  for(int i = 0; i < InFiles.dados.size();i++)
    Nsize+=InFiles.dados[i].size();


  // Input File Definition - IFD:  
  //     Number of fixed Parameters lines      -  IFD_NP
  //     Number of columns Parameters          - *IFD_P (array[IFD_NP])
  //     Spot Line(onde com. def. de spot - 1) -  IFD_SL
  //     Number of Spot Parameters             -  IFD_SP
  //     posicao da lei de obs. de Borda       -  IFD_B
  //     Numero de Colunas de Imagens          -  IFD_NCI

  const int IFD_NP = 3;
  const int IFD_P[IFD_NP] = {6, 7, 7};
  const int IFD_NCI = 12;
  //  const int IFD_SL = 3;
  //const int IFD_SP = 5;
  const int IFD_B = 2;

  DataFileReader fin(argv[2]);

  // vector<vector<double> > in_par(IFD_NP+1);


  InFiles.in_par.resize(IFD_NP+1);

  //Lendo parametros de entrada
  ++fin;
  for(int i = 0; i < IFD_NP; i++){
    if(fin.numberOfColumns() == IFD_P[i]){
      double data_in;
      for(int j = 0; j < IFD_P[i]; j++){
	 data_in = fin.readColumn<double>(j+1);
	if(i == 1 && j == 0) data_in *= M_PI/180;
	InFiles.in_par[i].push_back(data_in);
      }
    }
    else { 
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[2] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }
    //++fin;
    if(++fin) {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[2] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }
  }

  //Lendo coeficientes do obscurecimento de borda

  if(fin.numberOfColumns()==definicoes_obscurecimento(int(InFiles.in_par[2][IFD_B]))){
    double data_in;
    for(int i = 0; i < definicoes_obscurecimento(int(InFiles.in_par[2][IFD_B])); i++){
		data_in = fin.readColumn<double>(i+1);
      InFiles.in_par[3].push_back(data_in);
    }
  }
  else {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[2] << "\n\tLinha  : " << fin.lineNumber() << endl
	  << "\tExpecting " << definicoes_obscurecimento(int(InFiles.in_par[1][IFD_B])) <<endl
	  << "\tGot " << fin.numberOfColumns() << endl;
	  
    aborta(erro.str());
  }

  // Mais nada para ler, fechar arquivo...
  fin.close();

  // Nao esqueca de redefinir o nome de cada variavel!
  show_par(cout, InFiles.in_par);

  // Calcular ffac para o dado raio da estrela secundaria
	
	const int FFAC_COL = 5;
	const int FFAC_LINE = 2;
	double ffac = InFiles.in_par[0][5] > 0 ? Radii2Ffac(InFiles.in_par,FFAC_LINE,FFAC_COL) : InFiles.in_par[FFAC_LINE][FFAC_COL];
	InFiles.in_par[FFAC_LINE][FFAC_COL] = ffac;
	cout << InFiles.in_par[FFAC_LINE][FFAC_COL] << endl;
	//exit(0);
  // Defino classe aqui e passo vetor para estrutura InFiles
  // assim nao preciso redefinir Roche nos outros locais do
  // programa

  Roche estrela_local(int(InFiles.in_par[0][0]),org_par(InFiles.in_par));
  estrela_local.vorb(true);
  estrela_local.refaz_super();
	
  int _Np_ = ceil(InFiles.in_par[0][5] / ( InFiles.in_par[1][5]/1e5/(2.*M_PI*estrela_local.K2()) ) ); //Numero de perfis para combinar
	
	if (_Np_ < 2) { //This is as lower as it can go!
		_Np_ = 1;
	}
	
	cout << " Np = "<< _Np_ << endl;

	vector<double> vphase(InFiles.Fase.size()*_Np_);

	/*************
	 Note that the phases are organized such that the order of the pixel maps is retained. This way, pixel maps from
	 0 to InFiles.Fase.size() corresponds to the fases of the first cicle and so one. I am doing it like this so I
	 need not change the algorith for the derivatives for the chi-square.
	 *************/
	
	for (int nph = 0; nph < _Np_; ++nph) {
		for (int i_ph = 0; i_ph < InFiles.Fase.size(); ++i_ph) {
//			cout << i_ph+InFiles.Fase.size()*nph << "  " << i_ph << " " << _Np_ << " " << nph <<endl;
			vphase[i_ph+InFiles.Fase.size()*nph] =  _Np_ > 1 ? InFiles.Fase[i_ph] + InFiles.in_par[0][5]/(_Np_-1) * nph - InFiles.in_par[0][5]/2.0 : InFiles.Fase[i_ph];

		}
	}
//	exit(0);
  typedef vector<double>::const_iterator cit;
  int index=-1;
  vector<double> vimg(estrela_local.size());
  for(int i = 0; i < vimg.size(); i++) vimg[i] = 0.0;
	
  for(cit vf = vphase.begin(); vf != vphase.end(); vf++){

	  index = index >= InFiles.dados.size()-1 ? index = 0 : ++index;
	  
	  //cout << index << " " << *vf << " " << InFiles.dados.size() << endl;
	  
	vector<int> pixmap;
    estrela_local.pixmap(*vf,InFiles.dados[index],pixmap);
    InFiles.pixmap.push_back(pixmap);
	  
    for(int i = 0; i < vimg.size(); i++){
      if(pixmap[i] >= 0) vimg[i] = 1;
      //elif vimg[i] = 1.0;
    }
	  

  }

  ostringstream saida;

  saida << "# - Parametros orbitais:\n";

  for(int i = 0; i < InFiles.in_par.size(); i++){
    for(int j = 0; j < InFiles.in_par[i].size(); j++){
      if( i == 1 && j == 0) saida << InFiles.in_par[i][j]/M_PI*180 << ' ';
      else saida << InFiles.in_par[i][j] << " ";
    }
    saida << endl;
  }
	
  //1o. - Guardo o nome dos arquivos de entrada...
  //lis_in.rewind(); 

  DataFileReader lis_in2(argv[1]);

  saida << "# - Arquivos de entrada:\n";
  vector<double> sp_flux;
  while(!++lis_in2){   
    string input = lis_in2.readColumn<string>(1);
    saida << input << " ";
    //Calcular fluxo do espectro
    DataFileReader data(input);
    vector<vector<double> > dados(2);
    while(!++data){
      dados[0].push_back(data.readColumn<double>(1));
      dados[1].push_back(data.readColumn<double>(2));
    }
    //GSL_Interp<vector<vector<double> > > itr_data(dados,dados[0].size());
	  double nfac = 0;
	  for (int npt = 0; npt < dados[1].size(); ++npt) {
		  nfac+=dados[1][npt];
	  }
	  nfac/=dados[1].size();
    //sp_flux.push_back(itr_data.integrate());
	  sp_flux.push_back(nfac);
    data.close();
    //lis_in2.seteol();
    //++lis_in2;
  }
  saida << endl;
#ifndef OLD_VERSION
  saida << "# - Fluxo dos espectros:\n";
  for(int i = 0; i < sp_flux.size(); i++) saida<< sp_flux[i]<< " ";
#endif

  saida << endl
	<< "# - Mapeamento dos pixels:\n";

  for(int j = 0; j < InFiles.pixmap[0].size(); j++){
    saida << j;
    for(int i = 0; i < InFiles.pixmap.size(); i++)
      saida << " " << InFiles.pixmap[i].at(j);
    saida << endl;
  }
  saida.flush();

  ofstream file_out (argv[3]);
  file_out << saida.str();

  lis_in2.close();
  file_out.close();
  //escrevendo imagem de pixels viziveis I = 1 e pixels nao viziveis I = -1

  estrela_local.NewI(vimg);
  surface viz_img;
  int size = 50;
  estrela_local.openGL_Image(size, viz_img);

  ofstream imgs("lixo.map");

  if(!imgs){
    ostringstream erro;
    erro << "Arquivo " << "lixo.map" << " nao pode ser aberto!";
    aborta(erro.str());
  }

  //show_par(imgs, in_par);
  
  int Flag_pix = 0;
 
  for(int i = 0; i < viz_img.theta.size()-1; i+=2){
  //for (int i = 0; i < size; i++){
  // for (int j = 0; j < size*2;j+=2){
      imgs << viz_img.I[i] << " ";
      
    if(Flag_pix >= size) {
      Flag_pix = -1;
      imgs << endl;
      //cout << "-> ";
    }
    Flag_pix++;
    //cout << i << endl; //<<i*2*size+j << " " << i << " " << j <<endl;
  }

  imgs.close();

  return 0;
}


int definicoes_obscurecimento(const int law)
{
  vector<vector<double> > par;

  Light_Curve teste(par);
  
  return teste.npar_obs_borda(law);
}

void show_par(ostream& out, vector<vector<double> > &par)
{
  vector<vector<string> > par_def (4);

  const int p = 0;
  const int p0 = 1;
  const int p1 = 2;
  const int p2 = 3;

  par_def[p].push_back("SIZE = ");
  par_def[p].push_back("Xaim = ");
  par_def[p].push_back("FWHM = ");
  par_def[p].push_back("S/N = ");
  par_def[p].push_back("ABS_FLAG = ");
  par_def[p].push_back("texp/Porb = ");

  par_def[p0].push_back("i = ");
  par_def[p0].push_back("M1 = ");
  par_def[p0].push_back("M2 = ");
  par_def[p0].push_back("P = ");
  //par_def[p0].push_back("S/N = ");
  par_def[p0].push_back("Vmax = ");
  par_def[p0].push_back("lres = ");
  par_def[p0].push_back("Ires = ");
	
  par_def[p1].push_back("f2 = ");
  par_def[p1].push_back("grav_coef = ");
  par_def[p1].push_back("limb_law = ");
  par_def[p1].push_back("nspot = ");
  par_def[p1].push_back("T_pole = ");
  par_def[p1].push_back("ffac = ");
  par_def[p1].push_back("fd = ");

  for(int i = 0; i < par[p2].size(); i++){
    ostringstream str0;
    str0 << "x" << i << " = ";
    par_def[p2].push_back(str0.str());
  }

  out << "# Tiago Ribeiro - LinuxC++ - 2006\n#";

  for(int i = 0; i < par.size(); i++){
    int nc = 0;
    if(par[i].size() != par_def[i].size()){
      cout << "WARNING - Definicoes diferentes entre par e par_def.\n"
	   << "par = " << par[i].size() << endl
	   << "par_def = " << par_def[i].size() << endl;
      if(par[i].size() > par_def[i].size()) nc = par_def[i].size();
      else nc = par[i].size();
    }
    else nc = par[i].size();
    for(int j = 0; j < nc; j++){
      if( i < 4 ) out << par_def[i][j];
      else out << par_def[3][j];
      
      if( i == 1 && j == 0) out << par[i][j]*180/M_PI << ' ';
      else out << par[i][j] << ' ';
    }
    out << "\n#";
  }
  out << "V(Km/s)  Flux\n";
}

vector<vector<double> > org_par(vector<vector<double> > &in_par)
{

  vector<vector<double> > saida;

  // Infelismente nao existe maneira automatica de fazer isso. :(

  //  cout << endl;

  // Neste caso so preciso "deslocar" o vetor in_par de um e ja estou com o 
  // formato desejado!

  saida.push_back(in_par[1]); 
  saida.push_back(in_par[2]); 
  saida.push_back(in_par[3]); 
   
 return saida;
}

double Radii2Ffac(const vector<vector<double> > &in_par,int LIN,int COL)
{
	vector<vector<double> > par;
	
	for (int i = 1; i < in_par.size();i++){
		par.push_back(in_par[i]);
	}
	
	
	
	const double r2 = double(in_par[LIN][COL]);
	par[LIN-1][COL] = 1.0;
	cout <<r2 << " "<< par[LIN-1][COL] << endl;
	Roche estrela_local(3,par);
	// Fazendo minimizacao para achar ffac (procurando entre 1 e 10)
	
	double ffac2 = 10.;
	double ffac1 = 1.;
	double r2_mod1 = estrela_local.R2_rsun();
	estrela_local.ffac(ffac2);
	double r2_mod2 = estrela_local.R2_rsun();

	// Caso raio seja maior que o raio do lobo de roche faco o raio ser igual ao do lobo de roche
	if (r2_mod1 < r2){
		cout << "# - Selected Radius greater than maximum allowed radius... Using ffac = 1.0 (R2 = RL1 = " << r2_mod1 << " ) \n";
		return 1.0;
	}
	
	int niter = 0;
	int MAX_iter = 50;
	double r2_mod = r2_mod2;
	//cout << "# --> "<< fabs(r2_mod-r2)/r2 << endl;
	
	while ( ( ( fabs(r2_mod-r2)/r2 > 1e-3) && (niter < MAX_iter) ) ){
		
		double ffac = (ffac2+ffac1)/2.0;
		estrela_local.ffac(ffac);
		double r2_mod = estrela_local.R2_rsun();
		
		cout << "# --> "<< fabs(r2_mod-r2)/r2<< " " << ffac << " " << r2_mod << " " << r2_mod1 << " " << r2_mod2 <<" " << r2<<endl;
		
		if (fabs(r2_mod1-r2) > fabs(r2_mod2-r2)) {
			//cout << "Aqui 1\n";
			//estrela_local.ffac(ffac1);
			if (estrela_local.R2_rsun() >= r2){
			//	cout << "## - Should not have done this...\n" ;
				ffac1 = ffac;
				//estrela_local.ffac(ffac1);
				r2_mod1 = estrela_local.R2_rsun();				
				
			}
			else {
				ffac2 = ffac;	
				estrela_local.ffac(ffac);
				r2_mod2 = estrela_local.R2_rsun();
				
				
			}
			
			
		}
		else {
			//cout << "Aqui 2\n";
			//ffac2 = ffac;
			estrela_local.ffac(ffac);
			ffac2 = ffac;
			r2_mod2 = estrela_local.R2_rsun();


		}
		niter++;
		
	} 

	cout << "# [DONE]" << (ffac2+ffac1)/2.0 << " / " << niter << endl;
	
	return (ffac2+ffac1)/2.0;
	
}

